/**
 * @file Interpolation_2D.h
 * @author your name (you@domain.com)
 * @brief 
 * @version 0.1
 * @date 2020-10-19
 * 
 * @copyright Copyright (c) 2020
 * 
 */
#ifndef _INTERPOLATIONOPERATER_2D_H_
#define _INTERPOLATIONOPERATER_2D_H_

#include <iostream>
#include <cmath>
#include <vector>
#include "Miscellaneous.h"
namespace MyMultigrid_2D::Operator
{
    class InterpolationOperator_2D
    {
    public:
        std::vector<double> input;
        std::vector<double> output;
    public:
        InterpolationOperator_2D(){};
        ~InterpolationOperator_2D(){};
        virtual void DoInterpolation(int l) = 0;
        virtual void SetInput(std::vector<double> v) = 0;
        virtual std::vector<double> ReturnInterpolaionResult() = 0;
    };
    
    class LinearInterpolaionOperator_2D:public InterpolationOperator_2D
    {
    public:
        LinearInterpolaionOperator_2D(){};
        LinearInterpolaionOperator_2D(std::vector<double> v){LinearInterpolaionOperator_2D::input = v;};
        ~LinearInterpolaionOperator_2D(){};
        void DoInterpolation(int l)
        {
            int rhsize = (pow(2,l + 1) + 1);
            output.resize(rhsize * rhsize);
            double h_ = 1 / double(pow(2,l + 1));
            for(int k = 0;k < output.size();k++)
            {
                if(isboundary(k,l + 1) == true)
                {
                    int i = index_i(k,l + 1);
                    int j = index_j(k,l + 1);
                    if(j * h_ == 0 || j * h_ == 1)//上下边界；
                    {
                        if(j % 2 == 0)//粗网格存在该点；
                            output[k] = input[index(i / 2,j / 2,l)];
                        else
                            output[k] = 0.5 *(input[index(i / 2,j / 2,l)] + input[index(i / 2,j / 2 + 1,l)]);
                    }
                    else
                    {
                        if( i % 2 == 0)
                            output[k] = input[index(i / 2,j / 2, l)];
                        else
                            output[k] = 0.5 * (input[index(i / 2,j  /2 ,l)] + input[index(i / 2 + 1,j / 2,l)]);
                    }  
                }
            }
            for(int i = 0;i <= rhsize / 2 - 1;i++)
            {
                for(int j = 0;j <= rhsize / 2 - 1;j++)
                {
                    int v1 = index(2*i,2*j,l+1);
                    int v2 = index(2*i+1,2*j,l+1);
                    int v3 = index(2*i,2*j+1,l+1);
                    int v4 = index(2*i+1,2*j+1,l+1);
                    int u1 = index(i,j,l);
                    int u2 = index(i+1,j,l);
                    int u3 = index(i ,j + 1,l);
                    int u4 = index(i + 1,j + 1,l);
                    output[index(2*i,2*j,l+1)] = input[index(i,j,l)];
                    output[index(2*i+1,2*j,l+1)] = 0.5*(input[index(i,j,l)]+input[index(i+1,j,l)]);
                    output[index(2*i,2*j+1,l+1)] = 0.5*(input[index(i,j,l)]+input[index(i,j+1,l)]);
                    output[index(2*i+1,2*j+1,l+1)] = 0.25*(input[index(i,j,l)]+input[index(i+1,j,l)]+input[index(i,j+1,l)]+input[index(i+1,j+1,l)]);
                }
            }
        };
        void SetInput(std::vector<double> v){LinearInterpolaionOperator_2D::input = v;};

        std::vector<double> ReturnInterpolaionResult(){return output;};
    };  
}

#else

#endif
